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ABSTRACT 

We present detailed simulations of the interaction of a stream emanating from a mass- 
losing secondary with the core of a massive supergiant in the slow merger of the two 
stars inside a common envelope. The dynamics of the stream can be divided into 
a ballistic phase, starting at the L\ point, and a hydrodynamical phase where the 
stream interacts strongly with the core. Considering the merger of a 1 and 5M Q 
star with a 20 Mq evolved supergiant, we present two-dimensional hydrodynamical 
simulations using the PROMETHEUS code to demonstrate how the penetration depth 
and post-impact conditions depend on the initial properties of stream material (e.g. 
entropy, angular momentum, stream width) and the properties of the core (e.g. density 
structure and rotation rate). Using these results, we present a fitting formula for the 
entropy generated in the stream-core interaction and a recipe for the determination 
of the penetration depth based on a modified Bernoulli integral. 
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1 INTRODUCTION 

A large fraction of massive stars are known to be members of 
binary systems which are close enough to interact strongly 
via Roche-lobe overflow (see e.g. Garmany, Conti & Massey 
1980). One of the most dramatic, but poorly understood 
types of binary interactions involves the complete merger of 
the two components in a common envelope which leads to 
the formation of a single, rapidly rotating star with possi- 
bly some highly unusual properties. This type of evolution 
generally occurs as the result of dynamical mass transfer 
(Paczyhski & Sienkiewicz 1972) where a more massive star 
(usually a giant) overfills its Roche lobe by an ever increasing 
amount leading to mass transfer on a dynamical time-scale. 
Since the companion is only able to accrete matter on its 
thermal time-scale, most of the excess matter is believed 
to form a common envelope engulfing both the mass donor 
and the secondary. Inside this envelope, the core of the giant 
and the secondary form a close binary. Friction between this 
immersed binary and the envelope then causes the orbit to 
shrink. If the energy released in the process is sufficient to 
eject the envelope, this leads to the formation of a very close 
binary consisting of the core of the giant and a largely unper- 
turbed companion star (Paczyhski 1976; Taam & Sandqvist 
2000). However, if the envelope is not ejected in the early 
rapid spiral-in phase, this phase is followed by a phase of self- 
regulated spiral-in where all the frictional energy deposited 
in the envelope can be transported to the surface where it 



is being radiated away (Meyer & Meyer-Hofmeister 1979; 
Podsiadlowski 2001). For a massive primary with a mass of 
~ 20 Mo filling its Roche lobe as a red supergiant after the 
end of helium core burning, the characteristic time-scale for 
the self-regulated spiral-in phase is several 100 yr. However, 
at some point, the spiraling-in secondary itself will fill its 
critical tidal lobe and start to transfer mass to the core of 
the giant. This situation is in many respects similar to the 
case of normal mass transfer by Roche-lobe overflow except 
for several fundamental differences. (1) Mass transfer occurs 
inside a low-density, opaque envelope, where the density con- 
trast between the envelope and the secondary's average den- 
sity may be as high as 1:10 6 . (2) The mechanism for driving 
mass transfer is orbital angular momentum loss due to the 
friction of the immersed binary with the envelope. Because 
of the short spiral-in time-scale this leads to extremely high 
mass transfer rates of order 0.01 — lOMoyr . (3) Since 
the spiral-in time-scale is much shorter than a typical syn- 
chronization time-scale in a radiative star (e.g. Zahn 1975), 
the spin of the secondary may not be synchronized with the 
orbit. 

As a result of this interaction, the secondary is gradually 
being dissolved inside the common envelope. Even though 
the time-scale for this dissolution is short compared to the 
evolutionary time-scale of a giant, it is still much longer 
than the dynamical time-scale of the secondary. We there- 
fore refer to this process as a 'slow merger' of the compo- 
nents, although this phase is likely to be terminated by a 
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delayed dynamical instability when the flat-entropy core of 
the secondary is being exposed and mass transfer becomes 
catastrophic. This will ultimately lead to the dynamical dis- 
ruption of the remnant core (if it is relatively unevolved) . 

In this paper we present the first part of a systematic 
study of the slow merger of massive stars, concentrating on 
the interaction of the mass stream emanating from the mass- 
losing immersed star with the core of the giant. This stream 
will generally not intersect with itself and form an accretion 
disc, but will instead penetrate deep into the core of the 
giant before it is stopped by the increasing pressure inside 
the core. As we will show, the stream may be able to reach 
a radius as small as ~ 10 10 cm. Since this is deep inside the 
hydrogen-exhausted helium core of a 20 Mq star that has 
already exhausted helium in the core, this may have impor- 
tant consequences: (1) It can lead to the dredge-up of helium 
from the core which can then be mixed with the rest of the 
envelope. (2) Since the material in the stream is hydrogen- 
rich and is heated to temperatures as high as several 10 s K in 
the post-impact region, this may lead to some unusual nucle- 
osynthesis, in particular s-processing. How deep the stream 
can penetrate depends on the initial entropy of the stream 
material, but also its initial angular momentum and most 
importantly on the structure near the core and the amount 
of entropy generated in the interaction with the core. It is 
the purpose of this paper to systematically explore how the 
penetration depends on these various factors and to deter- 
mine the properties of the post-impact material. In a sub- 
sequent paper, we will use these results to model the slow 
merger of two massive stars in the context of the progen- 
itor of SN 1987A, whose highly anomalous properties are 
likely the result of such a merger (see Podsiadlowski 1992, 
1997 for discussion and references). While the hydrodynam- 
ical simulations presented in this paper were designed with 
this application in mind, our analysis is sufficiently general 
to allow application to other systems that have been sug- 
gested to be post-merger objects (e.g. V Hydrae [Kahane 
et al. 1996]; FK Comae stars [e.g. Rucinski 1990; Welty & 
Ramsey 1994]). 

In Section ^ we outline the basic assumptions for mod- 
elling the initial properties of the stream leaving the mass- 
losing secondary. In Section ^ we use this model to obtain a 
description of the stream at the start of the hydrodynami- 
cal simulations. Section ^ provides a description of the nu- 
merical method used in the hydrodynamical calculations. In 
Section |E| we present the results of detailed hydrodynami- 
cal simulations and systematically discuss the physics of the 
interaction between the jet and the ambient matter, in par- 
ticular the mechanism of the stream dissipation, the effects 
of rotation and of nuclear burning in the stream. Finally, in 
Section 6 we present an semi-empirical prescription which al- 
lows to estimate the penetration depth, taking into account 
the entropy generated in the interaction. 



2 THE INITIAL CONDITIONS OF THE 
STREAM 

2.1 Mass-loss rates 

Once the immersed secondary starts to fill its tidal lobe, 
it starts to transfer mass to the core of the giant. The 



mass-transfer rate is determined by the friction between the 
spiraling-in immersed binary and the envelope which drives 
the evolution of the system. Since the mass-loss time-scale 
is generally much shorter than the nuclear or thermal time- 
scale for a main-sequence star, the radius responds adiabat- 
ically to mass loss, i.e. 
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where R2 and Ma are the radius and the mass of the mass- 
losing secondary, respectively. 

The Roche-lobe radius J?rl of the secondary changes ac- 
cording to 
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where the Roche-lobe radius Rb.l is given by Eggleton 
(1983) 
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Here q — Ma /Mi is the mass ratio, and Mi is the total mass 
of the primary core within the binary separation D. 

Assuming that the radius of the star is equal to the 
Roche-lobe radius, we can equate equations (1) and (2) and 
obtain an estimate for the mass-loss rate from the secondary 
as 
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where d\nD/dt is the change in the orbital separation due 
to tidal and viscous friction. For typical conditions in a 
20 Mq supergiant, the resulting mass- loss rates are typi- 
cally in the range of 0.01 to lOMoyr - (Ivanova, Podsiad- 
lowski & Spruit 20001; and Ivanova & Podsiadlowski 2002 
[IP]). To obtain values for £ad> we performed separate stellar- 
evolution calculations where mass was taken off the assumed 
secondary on a time-scale much shorter than the thermal 
time-scale in its outer envelope (see e.g. Podsiadlowski, Rap- 
paport & Pfahl 2001). 



2.2 The stream density profile and the stream 
width 

In previous treatments of mass outflow from L\, it is gener- 
ally assumed that the stream is isothermal (see e.g. Lubow 
& Shu 1975; LS; Bisikalo et al. 1997). Because of the high 
expected mass- loss rate in our problem, it is more appropri- 
ate to assume that the stream behaves adiabatically, since 
the radiative diffusion time-scale is generally much longer 
than the characteristic flow time-scale (for a stream with a 
radius of 10 10 cm, the diffusion time is of order 10 6 s). 

In the case where the secondary's spin is aligned and 
synchronous with the orbit, we can determine the density 
distribution in the neighbourhood of the L\ point by ex- 
panding the potential around this point in a Taylor series 
(see e.g. LS), where we use Cartesian coordinates (x,y,z): 
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Here fi is the angular speed of the binary and fi — 
M 1 /(M 1 +M 2 ). 

A natural assumption is that the energy distribution 
within the stream, consisting of kinetic and internal ener- 
gies, follows the shape of the Roche potential and that the 
Bernoulli integral for the matter leaving the donor star is 
conserved (for an adiabatic stream) 
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Here ul 1 and cl x are the velocity and the sound speed of 
the matter at L\ , and 7 is the adiabatic index of the stream 
material. Let us assume that in the L\ neighbourhood the 
stream is laterally in hydrostatic equilibrium. This is a good 
approximation, especially in the Z— direction (LS). Then, in 
the Y Z— plane 
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where r is the distance to the Li point and a » A is a 
unitless constant which depends on the shape of the poten- 
tial in the neighbourhood of L\. For binary systems with a 
mass ratio, q, between 0.006 and 1.0, it has a value of 6.16 
to 8.0 (LS). Then the density distribution within an adia- 
batic stream can be obtained from the condition of lateral 
hydrostatic equilibrium as 
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where p c is the density of the material in the core of the 
stream, and W is the characteristic width of the stream, 
which depends on the properties (e.g. entropy, temperature) 
of the matter in the secondary: 

1/3 
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Here s = P/p~" is the adiabatic constant of the material from 
the secondary. An adiabatic stream can not have a width 
wider than W = W \j2/a (as can be seen from equ. ^|). 

We assume that the velocity of the material at L\ is 
given by the local sound speed (see LS) and that this speed 
varies through the stream cross-section. Then, integrating 
the mass-flow rate across the stream cross section near L\, 
we can relate the mass-loss rate M2 (in units of M© yr _1 ) to 
the density in the stream core for an adiabatic stream (with 
7 = 5/3): 



stant of the material from the secondary, ^4 = S7/10 4 s , 
x s = W/W is the unitless stream width. 

When the stream forms in the hot and high-pressure en- 
vironment of a common envelope, the actual stream width 
is determined by the local conditions of the envelope. If we 
assume that the outer edge of the stream is in pressure bal- 
ance with the common envelope, the stream density at the 
edge of the stream cannot be less than 
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where Pqe is the ambient pressure in the common envelope. 
The stream width for the current mass loss rate then can be 
found from the equ. (10) in combination with 

1/ / a ,\ -1/(7-1) 
Pc = (Pce/ S ) 1/7 (l - |^ 2 ) (12) 

(from equ. 8). 

2.3 The effect of non-synchronous rotation 

Since the secondary's spin may not be synchronous with the 
orbit (in fact, it may be rotating faster or more slowly de- 
pending on the circumstances), we also need to consider the 
case of non-synchronous rotation. For this purpose we adopt 
a coordinate system that rotates with the secondary. At the 
equipotential surface which determines the surface of the 
secondary (the donor), there is an area where matter is be- 
ing accelerated outwards from the stellar surface. From this 
zone matter may flow out. In the case of a synchronized sec- 
ondary, the physical situation is a close analogue of a Laval 
nozzle, where the equipotential surfaces take over the role 
of a nozzle. In the standard case, when a thin flow is con- 
sidered, Li can be taken as a point. Then, as for the case of 
the standard Laval nozzle, matter is almost completely unaf- 
fected by external acceleration at the nozzle mouth (at the 
sonic point), but experiences negative acceleration before 
and positive acceleration after passing though the nozzle. 

For an asynchronously rotating secondary, the position 
of the positive acceleration zone with respect to the sec- 
ondary changes with time. To parametrize the degree of a- 
synchronicity, we introduce a parameter /, 



/ = 



(13) 



where Q sec is the angular velocity of the secondary. This pa- 
rameter is zero if the secondary rotates synchronously. The 
position of the point where the acceleration experienced by 
matter on the stellar surface becomes non-negative can be 
found similarly to the L\ point (in order to avoid confusion 
with the standard L\ point, we refer to this point as the 
'leaving point', LP), considering the generalised Roche po- 
tential, for the frame with the origin at the center of the 
secondary, which rotates with the secondary (see e.g. Limber 
1963): 
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where S15 = (P/p 7 )/(10 15 cgs) is the scaled adiabatic con- 
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and (f, rj) are Cartesian coordinates in the frame rotating 
with the secondary. For a secondary which rotates faster 
than the binary itself, LP will normally be shifted from Li 
towards the donor; and for a secondary which rotates more 
slowly than the binary this point will be shifted towards the 
primary. A rough estimate of the linear speed with which 
the zone of positive acceleration travels across the secondary 
surface, ^zonc = fil sec r SBC , shows that it is typically twice the 
local sound speed of the secondary (for the case of / = 1 
and considering parameters characteristic for a system at the 
point where the secondary starts to overfill its Roche lobe). 
In this case, we do not have a normal nozzle (as for a syn- 
chronized secondary), as matter at a particular position on 
the secondary's surface (in the frame rotating with the sec- 
ondary) is not affected by any matter that has escaped from 
the secondary before it enters the zone of non-negative accel- 
eration. Therefore we may consider what happens to a fluid 
element on the secondary independently from any matter 
that has flowed out previously. A physical analogue of this 
situation is a ball, containing the (secondary) matter, that 
is separated from the outside by negative acceleration. At 
the surface of the ball there is a hole which travels with su- 
personic speed. As far as the matter in the ball is concerned, 
the surface of the ball acts like a wall that is suddenly be- 
ing removed (for some period of time) and then placed back 
again. Matter will escape from the ball through the hole 
only if the duration of the interval of positive acceleration 
is long enough for matter to be able to escape before it is 
pulled back by the negative acceleration. Another analogue 
of the problem is the sudden expansion of a gas sphere into 
a vacuum. In this case, the speed with which matter ex- 
pands can be as high as 3 times the sonic speed (see e.g. 
Zeldovich & Raizer 1966). In a real situation, we may ex- 
pect a combination of these and other effects that cannot 
be easily treated analytically. As this discussion shows, for 
an a-synchronously rotating secondary, the range of the pos- 
sible initial flow velocities includes significantly supersonic 
speeds. Fortunately, our ballistic calculations in § 3 show 
that the stream properties close to the primary core are not 
strongly dependent on the initial flow velocity in the stream. 

2.4 The escape angle - generalization of the LS 
formalism 

In the case where the secondary's spin is not synchronous 
with the orbit, the angle at which matter leaves the sec- 
ondary will be different from the synchronous case analyzed 
by LS. To obtain this angle in our problem, we generally fol- 
low the LS formalism with the following changes: we place 
the origin of the coordinate system not at L\ but at LP, 
and we use a coordinate frame which rotates with the sec- 
ondary. Using cylindrical coordinates, the potential in the 
neighbourhood of LP can be expanded (in the z = plane) 
as 
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where Xlp is the distance between the center of mass and 
LP and 8 is the angle between the flow direction of the 
stream and the direction of the primary. 



The escape angle of the stream in the frame rotating 
with the secondary can then be found, analogously to the 
derivation in LS, from the equation 
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For the case of a synchronously rotating secondary (f=0), 
this equation reduces to the solution for 8 S obtained by LS. 
The angle with which the stream will escape from LP in the 
frame co-rotating with the binary is given by the transfor- 
mation 
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where v e is the velocity of the stream at the nozzle (sonic 
for a synchronously rotating secondary). 



3 THE BALLISTIC PHASE 

The motion of the stream can be naturally divided into two 
phases, a (quasi-) ballistic phase where the stream trajectory 
is not strongly affected by the interaction with the ambient 
matter, which can be treated more or less ballistically, and a 
hydrodynamical phase where the dynamics of the stream is 
dominated by the hydrodynamical interaction with the core. 
Considering these two phases separately helps to avoid nu- 
merical problems (due to the complexities near the L\ point) 
at the start of the hydrodynamical calculation. However, it 
also introduces some uncertainty since matter in the inner 
region surrounding the primary core may not be in complete 
co-rotation with the orbit and will have a backward motion 
with respect to the co-rotating frame. Such relative rotation 
causes a force which reduces the interception angle between 
the falling stream and the normal direction to the core rela- 
tive to the predictions from a purely ballistic trajectory (see 
also § 5.3.1). 

The equations of motion of the stream are 
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Here u — \Ju% + Uy is the magnitude of the central stream 
velocity, and u x ,u y are its components. The pressure gradi- 
ent along the path (1/p) dP/dl is included for completeness. 
This term is generally small, since the Mach number of the 
flow is large except near the nozzle. It is kept here for com- 
pleteness and consistency in the calculations. The structure 
of the common envelope enters into these equations through 
the pressure boundary condition at the stream edge, which is 
set equal to the ambient pressure in the envelope. The den- 
sity within the stream is related to the pressure by the adia- 
batic constant of the secondary material s — P/p 1 = const. 
Together with a specified mass-flow rate, this set of equa- 
tions provides a complete set that completely describes the 
dynamics of the stream in this phase. 
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3.1 Initial conditions for the hydrodynamical 
calculations 

The initial conditions for the hydrodynamical simulations 
are taken directly from the parameters obtained from the 
ballistic calculations for a range of binary parameters. 
Specifically, we consider parameters representing the spiral- 
in and subsequent merger of a secondary of 1 and 5 Mq 
inside a 20 Mq evolved supergiant which has completed he- 
lium burning in the core and has a core mass of ~ 7Mq. 
To model the secondaries, we first followed their evolution 
to the same age as the primary and then evolved them fur- 
ther, subjecting them to very fast mass loss to model their 
adiabatic response to mass loss. These calculations give the 
radii of the stars and their surface entropies as a function of 
current stellar mass. 



4 THE NUMERICAL METHOD 
4.1 The hydrodynamics code 

For the hydrodynamical simulations of the stream-core in- 
teractions we use a code based on the PROMETHEUS hy- 
drodynamical code (Fryxell, Muller & Arnett 1989) . This is 
an Eulerian code, which uses a second-order Godunov-type 
scheme to solve the hydrodynamical equations, the piecewise 
parabolic method (PPM) of Colella & Woodward (1994). 
This code has been widely used by different groups (see e.g. 
Kercek, Hillebrandt & Truran 1998; KHT), and we refer to 
these papers for a detailed description of the code and nu- 
merical tests. To make it applicable to our problem, we had 
to make a number of modifications to the original code, in 
particular to the treatment of the gravity field, the boundary 
conditions, to the equation of state and the hydrodynamical 
equations themselves (to take into account the effects of a 
frame co-rotating with the binary). Since the mass transfer 
occurs in an opaque environment, radiative losses are not 
important and have been neglected. 

In the presence of a large pressure gradient at the core- 
envelope interface and a correspondingly strong gravita- 
tional field, Godunov-type schemes produce intrinsic accel- 
erations, which during a few dynamical time-scales td cre- 
ate significant outward motion (KHT). For a first-order Go- 
dunov method it is possible to find an analytical formula to 
modify the interface states. Then source terms (the grav- 
itational field in our case) will be balanced by flux dif- 
ferences (LeVeque 1998). This is not possible for higher- 
order Godunov-type schemes, where each problem requires 
a separate treatment (KHT). To reduce the artificial accel- 
eration in our case, we assume that the gravitational field 
can be considered as constant in time and that there is no 
self-gravity. Then, in each time step, we modify the inter- 
face values before applying the Riemann solver, by reducing 
them from both (left and right) interfaces by a pressure flux 
caused by gravity. Furthermore, we enforce the condition of 
hydrostatic equilibrium at the boundary (for ghost cellsFl) 



1 Ghost cells are a few additional cells (four in the case of a PPM 
code) at each boundary of the computational domain, where the 
values at these cells are determined by the boundary conditions. 



in the direction of the gravity field. These modifications en- 
sure that any initial model remains stable for arbitrarily long 
time. 



4.2 The model set-up 

We performed 2-dimensional calculations in polar cylindrical 
coordinates, where we usually used 300 x 300 grid points. In 
the case of very narrow streams, we increased the resolution 
to 600 x 600 grid points. We also carried out test calcula- 
tions with grids of higher resolution to satisfy ourselves that 
further increasing the number of grid points does not signif- 
icantly affect the results. As a general rule, we ensure that 
the inflowing stream contains at least 20 grid points in the 
azimuthal direction. We also carried out a few 3-dimensional 
calculations (in spherical coordinates) for the case of a sym- 
metrical (non-inclined) stream in a non-rotating frame. The 
results of these calculations showed that there were no sig- 
nificant differences between the 2- and 3-dimensional calcu- 
lations with the same initial parameters. In particular, the 
penetration depths of the stream were very similar in both 
cases. 



4-2.1 Ambient matter 

To model the region of interest for the stream-core in- 
teraction inside the common envelope, we adopt power- 
law distributions for the temperature and the pressure, i.e 
T(r) = T(r ) (r /r) aT and P(r) = P(r ) (r /r) ap . We fit- 
ted these power laws to the structure in actual CE calcu- 
lations (as described in detail in IP). Typical values for qt 
and ctp are in the range (1.3; 5.2) + (0.8; 3.2), but can be as 
high as (1.7; 7) (corresponding to a structure with an adia- 
batic index 7 am b ~ 1.44) and reasonably describe the regions 
of the stellar models at the the evolutionary stage of inter- 
est (IP). For our calculations here, we adopted parameters 
(ar; ctp) = (1.3; 5.2) and (0.8; 3.2) for models representing a 
20+1 and a 20+5 CE simulation, respectively. Throughout 
the domain of our calculations, we assume that the ambient 
matter has a constant composition, similar to the compo- 
sition in the core region of the primary, mainly helium. In 
the full stellar models, hydrogen is exhausted below a radius 
of 4 x 10 10 cm, and the hydrogen mass fraction increases to 
~ 0.1 at the outer edge of the domain (at 7.5 x 10 10 cm 
). We have tested that these differences do not affect the 
results appreciably. In this parametrized model, we assume 
that the ambient matter is initially in (quasi-)hydrostatic 
equilibrium; as the model is parametrized by the pressure 
and temperature gradient, the gravitational field is then de- 
fined implicitly by the initial pressure gradient of the frame 
and the initial density distribution. In calculations where nu- 
clear burning is included, we use a nuclear reactions network 
with 27 isotopes, which includes all important reactions up 
to 25 Mg, in particular all a reactions and the complete hot 
(/3-limited) CNO cycle. The reaction rates were taken from 
the Thielemann nuclear reactions library (provided by R. 
C. Cannon) The changes in chemical composition and the 
energy density due to the nuclear burning are updated after 
each time step in the hydrodynamical calculation. To pre- 
vent nuclear burning in the ambient medium, we assumed 
that the ambient medium contained neither H nor He in 
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these particular calculations. (For further details we refer to 
IP.) 



4-2.2 Boundary conditions 

The boundary conditions in the code, which describe the 
characteristics of the inflowing stream, are the internal and 
external Mach numbers, Mint and M ex t; the ratio between 
the central density in the stream and the ambient density 
at the outer boundary r\ p \ the initial inclination angle of the 
stream 9 (i.e. the angle between the initial flow direction 
and the direction of the primary). For the chemical compo- 
sition of the gas in the stream we use the abundances of the 
secondary. 

The boundary conditions for the rest of the box in the 
radial direction are outflow boundaries (except in the gas 
inflow zone) where the condition of hydrostatic equilibrium 
has been imposed for ghost cells. The boundary condition in 
the azimuthal direction also assumes outflow conditions (we 
used a grid which normally covers an angle ir). As Courant 
number we used A^fl = 0.6 -5- 0.8. A further parameter 
is the angular velocity with which the coordinate system 
rotates (most calculations were performed in the frame ro- 
tating with the primary core). 



5 INTERACTION BETWEEN THE STREAM 
AND THE AMBIENT MATTER 

5.1 The initial stream-core interaction: 
establishing a steady stream 

Figure ^ shows the initial interaction of a stream with a 
massive core and illustrates the development of a stationary 
stream. The parameters for this simulation are characteris- 
tic for a 1 Mq star filling its Roche lobe inside the envelope 
of a very evolved 20 Mq red supergiant (after helium core 
burning). The panels on the left show the hydrogen abun- 
dance, while the panels on the right show the divergence of 
the velocity field. The latter was chosen because it shows 
shock structures particularly clearly. 

The top panels show the stream just before its 'impact' 
with the core (569s after the beginning of the simulation), 
while the second set of panels is close to the point of impact 
(at 778s; note the bow shock in the right panels). Immedi- 
ately after the impact, stream matter bounces off the core, 
where most of it is just reflected by the core and continues to 
move in the forward direction (i.e. counter-clockwise), but 
some of it is pushed backwards, attaining an angular velocity 
component opposite to the rotation of the core. As material 
that has been stopped by the core is pushed by material in 
the stream following from behind, it starts to flow up again, 
driving two powerful shocks on the front and back side of 
the stream into the envelopes. These shocks compress the 
stream significantly. Once they have left the domain of the 
calculation, the stream has attained essentially a station- 
ary configuration, where the point of deepest penetration 
(R ~ 2 x 10 10 cm) and the stream shape no longer change 
significantly. The overall flow pattern also becomes more- 
or-less stationary, where all the matter leaving the stream 
in the core impact region flows up again, being vigorously 
mixed with helium from the core in the process. 



5.2 Mechanisms for stream dissipation 

The main objective of our calculation is to determine how 
deep the stream can penetrate into the core of a massive 
star rather than the initial transient behaviour. In the quasi- 
stationary situation this depends mainly on the entropy that 
is generated in the stream-core interaction. 



5.2.1 Pressure discontinuity and entropy generation 

While the stream is in lateral pressure balance with the am- 
bient matter at an early stage of infall, it becomes increas- 
ingly unbalanced as the velocity increases and becomes sig- 
nificantly supersonic relative to the ambient medium. The 
interaction between the stream and the ambient matter due 
to the jump in pressure can be treated in a simplified way 
as a Riemann problem (see e.g. LeVeque et al. 1998), a so- 
lution of which is a combination of shock and rarefaction 
waves. In the case where the stream expands, the rarefaction 
wave, propagating into the stream material, does not change 
the stream entropy. However, compression of the stream by 
the ambient matter, in the form of shocks moving into the 
stream, generates entropy inside the shocks. 

For a strong shock, where the pressure jump P am b /Pa 3> 
4, the coefficient of entropy change, i.e. the ratio of the en- 
tropy of the shocked stream material, S ss , to the initial en- 
tropy in the stream, S s , can be estimated as 



Ks = S ss /S s » const • P amb /P s 
(|Zeldovichl966|) . 



(21) 



In our case, we rather expect the development of weak 
shocks. Then an estimate for Ks can be written as 



Ks ~ 1 + const 



nT 1 ( 



Min 



(Pamb/Ps - 1) 



(22) 



( ^eldovichl966| ). The pressure jump Pamb/Ps changes with 
the distance to the primary core. In the case of a power-law 
parametrisation of the ambient pressure, equation (^2|) pre- 
dicts that models with smaller power-law indices will result 
in smaller pressure jumps for the same initial stream proper- 
ties and hence smaller values for Ks and that this coefficient 
increases as the power-law index increases. 

In Figure [| we present some of the key characteristics 
of a stream with the same parameters as in Figure [j], once 
a more-or-less stationary flow pattern has been established. 
Note in the top two panels that the entropy in the core 
of the stream is initially constant and that all the entropy 
generation takes places in two strong shocks at the stream 
boundaries (clearly seen in the V v panel). All of the stream 
matter has to pass through the intersection of these shock 
structures, where it essentially attains its final entropy and 
becomes comparable to the entropy of the ambient material, 
at least in the case of the log P/p 7 . For the entropy the differ- 
ence between ambient and stream matter is largely caused 
by the difference in molecular weights and its subsequent 
change is due to mixing of these two media rather than dissi- 
pation. At the point of deepest penetration (~ 2 x 10 10 cm), 
the stream matter has a temperature of 6 x 10 7 K and a 
density of ~ 3 gem -3 , high enough for nuclear burning (see 
§ 5.4). 
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Figure 1. The initial stream-core interaction and the establishment of a steady stream in a frame rotating with the core (with an 
angular velocity f2f ramc = 4 X 10~ 4 rad s — 1 ). The left panels show the distribution of hydrogen, the right panels the distribution of V • v. 
The initial Mach numbers in the stream, Mi nt , and ambient medium, M ex t are 8.5 and 3.2, respectively, the density ratio between the 
stream and the ambient medium, r) p = 8, the outer radius i?x = 7.5 X 10 10 cm, the angular radius of the stream R B = 0.037T and the 
initial stream inclination with respect to the radial direction = 13.5°. The velocity field is scaled according to log(t;/10 6 cms -1 ). 



8 N. Ivanova, Ph. Podsiadlowski and H. Spruit 




log Temperature lag Density 

Time 2039.7 3 Time Z039.7 s 

Figure 2. Characteristic properties (entropy, log P/p 1 , divergence of velocity, hydrogen mass fraction, density and temperature) for a 
stationary inclined stream (with the same parameters as in Fig. ^). The radius of deepest penetration is 2.0 X 10 10 cm in the simulation. 
At that depth the temperature in the stream is 6 X 10 7 K and the density 3gcm~ 3 . 
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Figure 3. Example of a Kelvin-Helmholtz instability experi- 
enced by a narrow stream with R s = 0.005-7T. (M- ln t = 8, M ex t = 
3.25, ri P = 8, R T = 7.5 X 10 10 cm, Qf rame = 4 X 10 _4 rad8 _1 , 
© = 18°.) The stream penetrates to a depth of 2.7 X 10 10 cm. 

5.2.2 Kelvin-Helmholtz instabilities and the stream width 

The K-H instability affects both the stream shape and the 
stream width. For small relative velocities, the stream is 



"sausage" -shaped and becomes increasingly "snake" -shaped 
as the relative velocity increases (see Fig. ^). It becomes 
more important as the stream narrows, since perturbations 
of a given wavelength are more important for a narrower 
flow. It therefore determines the minimum stream width and 
thus limits the depth of penetration. In particular, a stream 
that is very narrow initially will penetrate less deep than a 
stream with a similar parameters like velocity, density and 
pressure, but with a larger mass-flow rate and hence larger 
initial stream width. 



5.3 The effects of rotation 

5.3.1 A rotating medium 

In some cases, in particular near the end of the merging pro- 
cess, the envelope surrounding the spiraling-in components 
may no longer remain in co-rotation with the binary (ei- 
ther because the spiral-in time-scale becomes too short or 
the region surrounding the primary core expands rapidly). 
Then the stream feels an additional force due to the mov- 
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Figure 4. Examples for streams penetrating into a medium 
rotating relative to the binary. The stream rotates in an anti- 
clockwise and a clockwise direction relative to the medium in 
the top and bottom panels, respectively, with a relative velocity 
f^shift = 10 -4 rads - 1 . The initial positions of the stream centre in 
the azimuthal direction were taken to be 0.37T and 0.67T in the two 
cases, respectively. The other stream parameters are the same as 
in Figure | (M int = 8.5, M cxt = 3.2, n p = 8, R T = 7.5 X 10 10 cm, 
Rb = 0.03tt, Q amb = 4 X 10 -4 rads -1 , = 13.5°). The radii 
at the point of deepest penetration are 1.75 and 2.25 X 10 10 cm, 
respectively. 



ing external medium. This affects the stream's trajectory, in 
particular the angle of incidence, and thereby the penetra- 
tion depth. This is illustrated in Figure ^ which shows two 
stream calculations with the same parameters as in Figure [| 
but where the medium is assumed to rotate with an angular 
velocity fl B hut = 10 -4 rads -1 relative to the frame of the 
binary (rotating with angular velocity Q) either in the for- 
ward or in the backward direction. The frame of the medium 
(in which the calculation is performed) then rotates with a 
velocity S7 am b = S7±fi s hift. In this frame, the stream itself ro- 
tates either ant i- clockwise (top panel) or clockwise (bottom 
panel). As the stream is being pushed backwards by ambient 
matter moving against its direction of rotation, its trajectory 
steepens and consequently the penetration depth increases 
(top panel), while the opposite happens when the stream is 
pushed from behind (bottom panel). As a consequence the 
stream penetrates deeper to a radius of ~ 1.75 x 10 10 cm 
instead of ~ 2 x 10 10 cm in the standard case (Fig. ^), while 
it only reaches a depth of ~ 2.25 x 10 10 cm in the opposite 
case. The corresponding change in the stream temperature 
at the deepest point is only a few per cent. Note also that for 
the stream rotating anti-clockwise, the entropy generation at 
the back side of the stream, where the dynamical pressure is 
smaller than at the front side (and hence has a larger jump 
in pressure), is larger than for stationary inclined stream. 
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Figure 6. The logarithm of the adiabatic constant (s = P/p 7 , 
top panel) and specific angular momentum (bottom panel) near 
the centre of the streams as a function of penetration depth 
for the simulations in Figure |jj with angular velocities Q = 
(solid curves), 0.001 (dashed curves), 0.0025 (dot-dashed curves) 
and 0.005 rads -1 (dot-dot-dot-dashed curves). The correspond- 
ing angular-velocity profiles for the ambient core matter are 
shown in the bottom panel as lighter curves with the same line 
style. The thick solid curve in the top panel refers of logs in 
the ambient matter. The filled circles at the end of individual 
curves mark the points of deepest penetration. The ticks along 
each curve indicate particular values of the central stream hydro- 
gen abundance by mass (0.6, 0.5, 0.4, 0.3 from right to left). 



5.3.2 Core rotation 

As the stream interacts with the core, it injects not only 
matter but also angular momentum into the core, spinning 
it up in the process. Since the core only has to accrete a 
small fraction of the angular momentum available from the 
secondary to be spun up to critical rotation, we may expect 
that this itself will affect the dynamics of the interaction. To 
demonstrate this, we performed four simulations for streams 
with the same initial parameters and similar core properties 
except that we varied the rotation of the core from to 
0.005 rads -1 within a core radius of 4 x 10 10 cmQ. The results 
of these simulations are shown in Figures |and|[ 

As Figure [| shows, the core-stream interaction changes 
qualitatively and quantitatively, where the penetration 
depth decreases with increasing core rotation. In the case 
of a rapidly rotating core, the penetration depth is centrifu- 



2 Note that these cores are not completely comparable, since we 
assumed that the ambient core matter was in quasi-static dy- 
namical equilibrium and had the same pressure and temperature 
distributions in all cases. This implies slightly different gravita- 
tional accelerations (and hence different core masses). 
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Figure 5. Stream simulations illustrating the effects of rotation of the primary core (in a frame rotating with the binary). In the four 
calculations, the core is assumed to rotate with constant angular velocities of f2 cor c = (top left), 0.001 (top right), 0.0025 (bottom left) 
and 0.005 rads - 1 (bottom right), respectively, up to a radius 4 X 10 10 cm. The other stream parameters are M; nt = 9.8, M ext = 3.7, 
rip = 8.4, R T = 7.5 10 cm, ij s = 0.025tt, = 18°. 



gaily limited (i.e. by the angular momentum in the stream). 
In Figure |^ we show the logarithm of the adiabatic con- 
stant and the specific angular momentum near the center 
of the streams in the four simulations as a function of ra- 
dius. In the two simulations with no or slow core rotation 
(solid and dashed curves), the stream's angular velocity is 
efficiently being braked below a radius 2.6 x 10 10 cm by its 
interaction with the core and reaches the deepest point of 
penetration (marked as filled circles) where the adiabatic 
constant is comparable to the ambient one. The specific an- 
gular momentum is always substantially less than the local 
cr itical angular momentum for centrifugal support, given by 
ij gR 3 , where g is the gravitational acceleration at radius R. 
On the other hand, in the two cases with rapid core rotation 
(dot-dashed and dot-dot-dot-dashed curves), the stream is 
first accelerated in the azimuthal direction by the interac- 
tion with the core, in fact reaching an angular velocity that 
is larger than the local velocity of the core (even though the 
initial angular velocity at the core boundary was less) . At the 
point of deepest penetration, the specific angular momen- 
tum is close to the critical angular momentum, and hence 
the stream becomes centrifugally supported. Note also that, 
in this case, the stream does not bounce but starts to form 
a ring orbiting the core and gradually mixing with it. This 
demonstrates that the main effect of rapid core rotation is to 
prevent /reduce the azimuthal braking of the stream. While 
this leads to less entropy generation in the stream (see the 
top panel in Fig. |^), which ordinarily would mean that the 
stream should be able to penetrate deeper, it makes the an- 
gular momentum of the stream itself a key factor in limiting 
the penetration depth. In the simulation with the fastest 
core, the stream matter is still substantially denser than the 



ambient matter at the final point shown in the simulations. 
Since it is also rotating faster than the local core material, 
it is reasonable to expect that the stream will continue to 
spiral in as it is being braked by the ambient medium and 
will penetrate somewhat deeper than is shown in the figures. 
(We were unable to follow the subsequent evolution since the 
stream matter was leaving the domain of the calculation at 
this point.) 

Even in the case where the core is non-rotating, the 
specific angular momentum in the stream increases during 
a portion of its trajectory (see the solid curve in the bot- 
tom panel between 3 and 3.8 x 10 10 cm) . This is a direct 
consequence of the stream splitting into two components in 
the impact region, a forward and a backward component 
(this can be seen more clearly in Figs 1 and 2). As mat- 
ter in the forward component is pushed by matter following 
from behind, it is being accelerated and gains angular mo- 
mentum, while matter flowing backwards attains negative 
angular momentum. 

Finally, note that because of the narrower stream in the 
post-impact region of the more rapidly rotating core, stream 
material mixes more efficiently with the core material (Fig. || 
and ^). This may have been important consequences for the 
nucleosynthesis in this region (see § 5.4 and IP). 

5.4 The effects of nuclear burning 

One of the objectives of our study is to investigate the nu- 
cleosynthesis in the stream material after it has been heated 
by its impact with the core of the primary. The associated 
heating by nuclear burning itself may then affect the hydro- 
dynamics of the stream. As can already be seen from the 
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Figure 7. Density and temperature distributions for stream simulations illustrating the effect of violent nuclear ignition. The top two 
panels show the stream just before impact with the core and the point of nuclear ignition. The middle two panels show the detonation 
wave that develops shortly after nuclear ignition that disrupts the inflowing stream. The bottom two panels show a stream at the same 
time if nuclear burning is not included for comparison. The stream parameters are M; nt = 9.9, M cx t = 4., r\ p = 7, i?x = 5 X 10 10 cm, 
Rs = 0.03tt, 6 = 10°. 



simulations presented so far, the characteristic time-scale 
for the stream infall is of order 1000 s, where the stream 
material typically spends no more than 50 — 100 s near the 
impact region before it bounces or forms a ring orbiting the 
core (in the case of a fast rotating core). Nuclear burning can 
only affect the stream motion if the characteristic time-scale 
for nuclear burning is of order or shorter than the stream 
dynamical time-scale. For the largest part of the stream tra- 
jectory, the stream material is relatively cold (for the range 
of stream models considered). 

Nuclear burning only starts to increase dramatically af- 
ter it has been heated to a temperature comparable to the 
ambient temperature, which occurs suddenly in the region 
of impact. However, since the characteristic nuclear time- 
scale is generally longer than the local dynamical time-scale, 
stream material will have already moved far away from the 
region of impact and hence have little effect on the stream 
dynamics (in the pre-impact region) and the stream trajec- 
tory. It mainly affects the properties of the ambient mat- 



ter (its composition and temperature) and hence the back- 
ground conditions for the stream on a time-scale that is long 
compared to the stream infall time-scale (these effects and 
the associated nucleosynthesis can be followed with a stan- 
dard stellar evolution code; see IPS for a description of the 
nucleosynthesis calculations) . 

The situation can be different when the stream im- 
pacts the core more-or-less radially and when the core has a 
very steep density gradient, causing a very strong shock and 
hence very high post-shock temperatures. This can cause a 
nuclear runaway as illustrated in Figure (?]. As the middle 
panels in this figure show, it triggers a detonation wave that 
propagates up into the stream, almost completely disrupt- 
ing it in the process. The expansion wave propagating into 
the ambient medium causes a dramatic expansion of the en- 
velope material. Once all of the matter burning in this run- 
away has been consumed, one may expect that the stream 
re-establishes itself and that the process repeats, leading to 
pulsed nuclear burning. 
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6 THE PENETRATION DEPTH 
6.1 The entropy change 

In the idealized case of a stream of constant entropy, the 
depth of penetration can easily be estimated from the bal- 
ancing of the stream ram pressure with the ambient pres- 
sure. In our case, the pressure of the stream material be- 
comes comparable to the ambient pressure well before the 
stream reaches its deepest point. In this case it is better to 
use an estimate based on the conservation of the Bernoulli 
integral for matter along a streamline (though both esti- 
mates would give similar results) . However, since the stream 
entropy changes because of the compression by shock waves 
propagating into the stream from the sides, one has to use 
a modified Bernoulli integral to estimate the penetration 
depth. This is possible as long as the stream is wide enough 
that it is not strongly affected by the K-H instability and 
that core rotation plays a minor role. 

For this purpose, we rewrite the Bernoulli integral for 
stream material in the form: 
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(23) 



where u is the velocity of stream material, c the sonic speed 
of the stream material and 4> the potential. The subscript S 
corresponds to matter in the stream, the subscript T refers 
to the initial conditions and the subscript B to the con- 
ditions at the point of deepest penetration. The point of 
deepest penetration can then be determined by the condi- 
tion that the pressure of the stream is equal to the pressure 
in the ambient matter and that the final stream velocity 
is zero (consistent with the results of our hydrodynamical 
calculations) . 

For the stellar structure parametrized by the pressure 
gradient the corresponding potential can be written as 



$(r) = -Q P 
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(24) 
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where the subscript A corresponds to the ambient matter at 
the initial point of the stream. Assuming that the entropy 
of the stream changes according to 
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one can rewrite the Bernoulli integral to find K s as a func- 
tion of relative penetration depth, 8 = r/r-r, and the initial 
conditions for the stream and ambient matter: 
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Figure 8. The entropy change as the function of the Mach num- 
bers ratio for a parametrized ambient structure with ap = 5.2 
and a-p = 1.3. Asterisks correspond to values of Kg obtained 
from the hydrodynamical simulations by determining the appro- 
priate value for K B that would give the same penetration depth 
when using the Bernoulli integral as in the actual hydrodynami- 
cal simulation. The error bars are obtained by assuming an error 
in the determination of the penetration depth of 3 per cent. Solid 
curves correspond to the theoretical fit with k = 0.45, while the 
dashed curves show the change in entropy of the stream material 
in the hydrodynamical simulations. 



structures of the ambient medium, (qp;ot) = (5.2.1.3) and 
(4.2; 0.9), representing typical profiles in massive CE sys- 
tems. We then calculated the stream penetration depth for 
3 values of density ratios, rj p = 2, 10 and 30, and for a num- 
ber of models with different Mach number ratios, Mi n t/M GX t 
(where we varied Mint, keeping M cx t fixed). All streams were 
calculated with an initial stream inclination of 18° . 

Using the expression for a weak shock (equ. 25), we fit 
the results to 



1 + fc(ap, qt) ■ rj 



r ( 



Mi, 



(27) 



For an ambient medium, characterized by ap = 5.2; «t = 
1.3, these results are well fitted by a constant k = 0.45 and 
for a structure with ap = 4.2; «t = 0.9 by k = 0.1 (see 
Fig.|). 
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Here pA is the density of the ambient material at the initial 
radius of the steam, and /3 is the ratio of gas pressure to 
total pressure for the ambient matter. 

6.2 Numerical fits 

To determine the entropy change factor, K B , we performed 
sets of hydrodynamical simulations for two parametrized 



7 CONCLUSIONS 

We have investigated the physics of the stream interaction 
with the core of a massive star in the phase where a star 
immersed in a common envelope merges with the core and 
systematically explored how the penetration depth depends 
on the initial properties of the stream and the core struc- 
ture. For simulations representing the merger of a 1 and 
5Mq with an evolved 20 Mq supergiant, it typically pene- 
trates to a depth of 1 — 2 x 10 10 cm, where it is being mixed 
vigorously with core material and is heated to temperatures 
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> 10 s K. In the case of a slowly rotating core, the penetra- 
tion depth mainly depends on the initial entropy and the 
entropy that is generated in the interaction. This also de- 
pends on the initial width of the stream (which is mainly a 
function of the mass-transfer rate) , since narrow streams are 
unstable to a Kelvin-Helmholtz instability. The main role of 
angular momentum in the stream is to determine the angle 
of incidence of the stream (which affects the dissipation in 
the stream). However, for a slowly rotating core, it is effi- 
ciently braked by the interaction with the core, spinning up 
the core in the process. In the case of a rapidly rotating core, 
the stream gains angular momentum from the core and the 
penetration is centrifugally limited. The stream behaviour is 
qualitatively and quantitatively different in the two cases. If 
the core is rotating slowly, the stream bounces off the core, 
quickly expanding away from the core, while in the case of 
a rapidly rotating core, the stream merges with the core in 
the form of a centrifugally supported ring. 

In most cases, nuclear burning in the stream, which only 
becomes effective near the region of impact, does not affect 
the structure and dynamics of the stream itself, but affects 
the ambient medium which affects the long-term evolution 
of stream-core interaction. On the other hand, in the case 
of a very strong shock (in the case of a steep density gradi- 
ent in the core and more or less normal incidence), shock- 
triggered nucleosynthesis can trigger a detonation wave that 
propagates upstream, temporarily disrupting the stream. 

Using a simple model for the dissipation in a weak 
shock, we derived a simple recipe that allows the approx- 
imate determination of the penetration depth for a wide 
range of conditions. In a follow-up paper (IP), we will im- 
plement these results in a modified stellar-evolution code 
to model the whole merging phase of two stars in a slow 
merger. 
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